Integrando SAGA con R para el análisis de terreno
Antony Barja
2020-01-19
Introducción
Es fascinante poder conocer ciertos aspectos del mundo del software libre y de código abierto orientado al análisis de datos espaciales, en está oportunidad quiero mostrarte un pequeño entrenamiento integrando dos softwares super potentes dentro del campo de la geocomputación, el primero es SAGA, el cual es un software de Sistema de Información Geográfica de código abierto y el otro es R, el cual es un lenguaje de programación libre e interpretado y orientado a objetos, este mismo presenta un gran ecosistema de paquetes para una diversidad de ejes temáticos, entre ellas podemos mencionar a la estadística, la estadística espacial, la teledetección, web mapping, machine learning, entre otros.
Para este entrenamiento tenemos que tener ciertos requerimientos como tener instalado SAGA y R, de preferencia se recomienda usar la IDE de Rstudio para tener una forma más amigable al momento de programar con R.
Observación:
Se recomienda tener instalado la versión de SAGA 2.3 LTS hasta la 7.0 .
Si tienes instalado QGIS no es necesario instalar SAGA ya que viene por defecto dentro su ecosistema.
- Es importante tener en claro dos conceptos básicos al momento de trabajar con SAGA, estos son las liberías y módulos, entonces:
- Liberías: Hace referencia a todas los “ejes temáticos” que puedes trabajar en SAGA.
- Módulos: Hace referencia a todo los algoritmos de estos ejes temáticos.
Por último, para integrar SAGA con R tenemos que instalar el paquete RSAGA (Alexander Brenning et.al,2018), posteriormente debemos identificar el directorio de trabajo donde está alojado la ruta de los módulos de SAGA.
Instalación de los paquetes necesarios para nuestro entrenamiento
libs <- c('RSAGA', 'tidyverse','raster',
'sp','sf', 'tmap', 'cptcity')
install.packages(libs, dependencies = TRUE)Activación de los paquetes instalados
library(tidyverse) # Pkgs para ciencia de datos
library(RSAGA) # Pkgs para integrar SAGA con R
library(raster) # Pkgs para manejar datos raster
library(raster) # Pkgs para manejar datos vectoriales
library(tmap) # Pkgs para elaborar mapas temáticos
library(cptcity) # Pkgs para obtener paleta de coloresIndentificando el directorio de trabajo donde están los módulos de SAGA
Observación:
- Si usas un sistema operativo Windows y tienes instalado QGIS solo tienes que identificar la siguiente ruta:
env <- rsaga.env(path = 'C:/Program Files/QGIS 3.10/apps/saga-ltr')- Si usas el sistema operativo de GNU/Linux sólo necesitas usar la siguiente expresión:
env <- rsaga.env()## Search for SAGA command line program and modules...
## Doneenv## $workspace
## [1] "."
##
## $cmd
## [1] "saga_cmd"
##
## $path
## [1] "/usr/bin"
##
## $modules
## [1] "/usr/lib/x86_64-linux-gnu/saga"
##
## $version
## [1] "2.3.1"
##
## $cores
## [1] NA
##
## $parallel
## [1] FALSE
##
## $lib.prefix
## [1] "lib"El proceso anterior nos ayuda a integrar SAGA con R, adicionalmente nos brinda una información acerca de la versión de SAGA con la cual estamos trabajando, la ruta de los módulos y alguna información adicional de interés.
Ahora vamos a proceder a visualizar las librerias disponibles de SAGA para trabajar dentro de R
libs <- rsaga.get.libraries(path = env$modules)
libs## [1] "climate_tools" "contrib_perego"
## [3] "db_odbc" "db_pgsql"
## [5] "docs_html" "docs_pdf"
## [7] "garden_3d_viewer" "garden_fractals"
## [9] "garden_games" "garden_learn_to_program"
## [11] "garden_webservices" "grid_analysis"
## [13] "grid_calculus_bsl" "grid_calculus"
## [15] "grid_filter" "grid_gridding"
## [17] "grid_spline" "grid_tools"
## [19] "grid_visualisation" "imagery_classification"
## [21] "imagery_maxent" "imagery_opencv"
## [23] "imagery_photogrammetry" "imagery_segmentation"
## [25] "imagery_svm" "imagery_tools"
## [27] "imagery_vigra" "io_esri_e00"
## [29] "io_gdal" "io_gps"
## [31] "io_grid_image" "io_grid"
## [33] "io_shapes_dxf" "io_shapes"
## [35] "io_table" "io_virtual"
## [37] "pj_georeference" "pj_proj4"
## [39] "pointcloud_tools" "pointcloud_viewer"
## [41] "shapes_grid" "shapes_lines"
## [43] "shapes_points" "shapes_polygons"
## [45] "shapes_tools" "shapes_transect"
## [47] "sim_cellular_automata" "sim_ecosystems_hugget"
## [49] "sim_erosion" "sim_hydrology"
## [51] "sim_ihacres" "sim_qm_of_esp"
## [53] "sim_rivflow" "statistics_grid"
## [55] "statistics_kriging" "statistics_points"
## [57] "statistics_regression" "ta_channels"
## [59] "ta_compound" "ta_hydrology"
## [61] "ta_lighting" "ta_morphometry"
## [63] "ta_preprocessor" "ta_profiles"
## [65] "ta_slope_stability" "table_calculus"
## [67] "table_tools" "tin_tools"
## [69] "tin_viewer"# Usamos la función lenght para indentificar el total de liberías:
length(libs)## [1] 69Ahora vamos extraer todos los módulos de las liberías usando un pequeño código gracias al paquete tidyverse
# Estos son los módulos para una libería específica ('ta_channels')
mod <- rsaga.get.modules(libs = 'ta_channels', env = env)
# Ahora vamos a crear una tabla con toda las liberías
# y sus respectivos módulos:
rsaga_mod <- rsaga.get.modules(libs = libs,env = env) %>%
bind_rows(.id = 'librerias') %>% as_tibble()
# Solo visualizamos las 10 primeras filas de toda la tabla
head(rsaga_mod,10)## # A tibble: 10 x 4
## librerias code name interactive
## <chr> <dbl> <chr> <lgl>
## 1 climate_tools 0 Multi Level to Surface Interpolation FALSE
## 2 climate_tools 1 Multi Level to Points Interpolation FALSE
## 3 climate_tools 2 Earth's Orbital Parameters FALSE
## 4 climate_tools 3 Annual Course of Daily Insolation FALSE
## 5 climate_tools 4 Daily Insolation over Latitude FALSE
## 6 climate_tools 5 Monthly Global by Latitude FALSE
## 7 climate_tools 6 PET (after Hargreaves, Table) FALSE
## 8 climate_tools 7 Daily to Hourly PET FALSE
## 9 climate_tools 8 PET (after Hargreaves, Grid) FALSE
## 10 climate_tools 9 Sunrise and Sunset FALSEAhora vamos a trabajar con el paquete raster para poder leer el Modelo de Elevación Digital (DEM) que está esta disponible en el siguiente enlace click aquí
dem <- raster('Insumos/DEM.tif')
plot(dem)
Es importante corregir algunos valores vacíos que presenta nuestro modelo de elevación digital, para esto vamos a usar el algoritmo de Fill sinks de Wang & Liu.
rsaga.fill.sinks(in.dem = 'Insumos/DEM.tif',
out.dem = 'Insumos/DEM_fill')## Search for SAGA command line program and modules...
## Done
##
##
## SAGA Version: 2.3.1
##
## library path: /usr/lib/x86_64-linux-gnu/saga/
## library name: libta_preprocessor
## library : Preprocessing
## tool : Fill Sinks (Planchon/Darboux, 2001)
## author : Copyrights (c) 2003 by Volker Wichmann
## processors : 4 [4]
##
## loading: DEM
##
##
## Parameters
##
##
## Grid system: 30; 558x 661y; 963390.8124x 8058312.3304y
## DEM: DEM
## Filled DEM: Filled DEM
## Minimum Slope [Degree]: 0.010000
##
## Save grid: Insumos/DEM_fill.sgrd...# Visualización de nuestro DEM corregido
tmap_mode('view')
raster('Insumos/DEM_fill.sdat') %>%
tm_shape() +
tm_basemap('Esri.WorldTopoMap') +
tm_raster(style = 'cont',alpha = 0.75,
palette = cpt(pal = 'ncl_topo_15lev'))Ahora procedemos a calcular el mapa de pendientes, aspecto y de curvatura
rsaga.slope.asp.curv(in.dem = 'Insumos/DEM_fill.sdat',
out.slope = 'Salida/pendiente',
out.aspect = 'Salida/aspecto',
out.cgene = 'Salida/curvatura',
unit.slope = 1,
unit.aspect = 1,
method = 'maxslope')## Search for SAGA command line program and modules...
## Done
##
##
## SAGA Version: 2.3.1
##
## library path: /usr/lib/x86_64-linux-gnu/saga/
## library name: libta_morphometry
## library : Morphometry
## tool : Slope, Aspect, Curvature
## author : O.Conrad (c) 2001
## processors : 4 [4]
##
## loading: DEM_fill
##
##
## Parameters
##
##
## Grid system: 30; 558x 661y; 963390.8124x 8058312.3304y
## Elevation: DEM_fill
## Slope: Slope
## Aspect: Aspect
## General Curvature: General Curvature
## Profile Curvature: <not set>
## Plan Curvature: <not set>
## Flow Line Curvature: <not set>
## Method: maximum slope (Travis et al. 1975)
## Slope Units: degree
## Aspect Units: degree
##
## Save grid: Salida/pendiente.sgrd...
## Save grid: Salida/aspecto.sgrd...
## Save grid: Salida/curvatura.sgrd...# Ahora procedemos a visualizar el mapa de pendiente
raster('Salida/pendiente.sdat') %>%
tm_shape() +
tm_basemap('Esri.WorldTopoMap') +
tm_raster(style = 'quantile',n = 7,palette = 'Spectral')# Visualizando el mapa de aspectos
raster('Salida/aspecto.sdat') %>%
tm_shape() +
tm_basemap('Esri.WorldTopoMap') +
tm_raster(style = 'cont',palette = cpt(pal = 'grass_gdd',rev = TRUE))# Visualizando el mapa de curvatura general
raster('Salida/curvatura.sdat') %>%
tm_shape() +
tm_basemap('Esri.WorldTopoMap') +
tm_raster(style = 'cont',palette = cpt(pal = 'grass_ramp'))Ahora vamos a calcular el índice de humedad topográfica (TWI) y el índice de posición topográfica (TPI)
rsaga.get.usage(lib = 'ta_hydrology',module = 15)## Search for SAGA command line program and modules...
## Done
## library path: /usr/lib/x86_64-linux-gnu/saga/
## library name: libta_hydrology
## library : Hydrology
## Usage: saga_cmd ta_hydrology 15 [-DEM <str>] [-WEIGHT <str>] [-AREA <str>] [-SLOPE <str>] [-AREA_MOD <str>] [-TWI <str>] [-SUCTION <double>] [-AREA_TYPE <str>] [-SLOPE_TYPE <str>] [-SLOPE_MIN <double>] [-SLOPE_OFF <double>] [-SLOPE_WEIGHT <double>]
## -DEM:<str> Elevation
## Grid (input)
## -WEIGHT:<str> Weights
## Grid (optional input)
## -AREA:<str> Catchment area
## Grid (output)
## -SLOPE:<str> Catchment slope
## Grid (output)
## -AREA_MOD:<str> Modified Catchment Area
## Grid (output)
## -TWI:<str> Topographic Wetness Index
## Grid (output)
## -SUCTION:<double> Suction
## Floating point
## Minimum: 0.000000
## Default: 10.000000
## -AREA_TYPE:<str> Type of Area
## Choice
## Available Choices:
## [0] absolute catchment area
## [1] square root of catchment area
## [2] specific catchment area
## Default: 1
## -SLOPE_TYPE:<str> Type of Slope
## Choice
## Available Choices:
## [0] local slope
## [1] catchment slope
## Default: 1
## -SLOPE_MIN:<double> Minimum Slope
## Floating point
## Minimum: 0.000000
## Default: 0.000000
## -SLOPE_OFF:<double> Offset Slope
## Floating point
## Minimum: 0.000000
## Default: 0.100000
## -SLOPE_WEIGHT:<double> Slope Weighting
## Floating point
## Minimum: 0.000000rsaga.geoprocessor(lib = 'ta_hydrology',module = 15,
param = list(DEM = 'Insumos/DEM_fill.sdat',
TWI = 'Salida/Twi' ),
env = env)##
##
## SAGA Version: 2.3.1
##
## library path: /usr/lib/x86_64-linux-gnu/saga/
## library name: libta_hydrology
## library : Hydrology
## tool : SAGA Wetness Index
## author : (c) 2001 by J.Boehner, O.Conrad
## processors : 4 [4]
##
## loading: DEM_fill
##
##
## Parameters
##
##
## Grid system: 30; 558x 661y; 963390.8124x 8058312.3304y
## Elevation: DEM_fill
## Weights: <not set>
## Catchment area: Catchment area
## Catchment slope: Catchment slope
## Modified Catchment Area: Modified Catchment Area
## Topographic Wetness Index: Topographic Wetness Index
## Suction: 10.000000
## Type of Area: square root of catchment area
## Type of Slope: catchment slope
## Minimum Slope: 0.000000
## Offset Slope: 0.100000
## Slope Weighting: 1.000000
##
## Create index: DEM_fill
## catchment area and slope...
## pass 1 (39594 > 0)
## pass 2 (9653 > 0)
## pass 3 (4059 > 0)
## pass 4 (1920 > 0)
## pass 5 (950 > 0)
## pass 6 (484 > 0)
## pass 7 (246 > 0)
## pass 8 (130 > 0)
## pass 9 (64 > 0)
## pass 10 (32 > 0)
## pass 11 (18 > 0)
## pass 12 (5 > 0)
## pass 13 (2 > 0)
## pass 14 (0 > 0)
## post-processing...
## topographic wetness index...
## Save grid: Salida/Twi...rsaga.get.usage(lib = 'ta_morphometry',module = 18)## Search for SAGA command line program and modules...
## Done
## library path: /usr/lib/x86_64-linux-gnu/saga/
## library name: libta_morphometry
## library : Morphometry
## Usage: saga_cmd ta_morphometry 18 [-DEM <str>] [-TPI <str>] [-STANDARD <str>] [-RADIUS_MIN <double>] [-RADIUS_MAX <double>] [-DW_WEIGHTING <str>] [-DW_IDW_POWER <double>] [-DW_IDW_OFFSET <str>] [-DW_BANDWIDTH <double>]
## -DEM:<str> Elevation
## Grid (input)
## -TPI:<str> Topographic Position Index
## Grid (output)
## -STANDARD:<str> Standardize
## Boolean
## Default: 0
## -RADIUS_MIN:<double> Radius
## Value range
## -RADIUS_MAX:<double> Radius
## Value range
## -DW_WEIGHTING:<str> Weighting Function
## Choice
## Available Choices:
## [0] no distance weighting
## [1] inverse distance to a power
## [2] exponential
## [3] gaussian weighting
## Default: 0
## -DW_IDW_POWER:<double> Inverse Distance Weighting Power
## Floating point
## Minimum: 0.000000
## Default: 1.000000
## -DW_IDW_OFFSET:<str> Inverse Distance Offset
## Boolean
## Default: 1
## -DW_BANDWIDTH:<double> Gaussian and Exponential Weighting Bandwidth
## Floating point
## Minimum: 0.000000rsaga.geoprocessor(lib = 'ta_morphometry',module = 18,
param = list(DEM = 'Insumos/DEM_fill.sdat',
TPI = 'Salida/Tpi' ),
env = env)##
##
## SAGA Version: 2.3.1
##
## library path: /usr/lib/x86_64-linux-gnu/saga/
## library name: libta_morphometry
## library : Morphometry
## tool : Topographic Position Index (TPI)
## author : O.Conrad (c) 2011
## processors : 4 [4]
##
## loading: DEM_fill
##
##
## Parameters
##
##
## Grid system: 30; 558x 661y; 963390.8124x 8058312.3304y
## Elevation: DEM_fill
## Topographic Position Index: Topographic Position Index
## Standardize: no
## Radius: 0.000000; 100.000000
## Weighting Function: no distance weighting
##
## Save grid: Salida/Tpi...# Visualizando el mapa de índice de humedad topográfica
raster('Salida/Twi.sdat') %>%
tm_shape() +
tm_basemap('Esri.WorldTopoMap') +
tm_raster(style = 'cont',palette = cpt(pal = 'idv_blues'))# Visualizando el mapa de índice de posición topográfica
raster('Salida/Tpi.sdat') %>%
tm_shape() +
tm_basemap('Esri.WorldTopoMap') +
tm_raster(style = 'cont',palette = cpt(pal = 'neota_base_yellow'))** Con RSAGA tambien podemos el orden de Strahler y la red de drenaje**
rsaga.get.usage(lib = 'ta_channels',module = 5)## Search for SAGA command line program and modules...
## Done
## library path: /usr/lib/x86_64-linux-gnu/saga/
## library name: libta_channels
## library : Channels
## Usage: saga_cmd ta_channels 5 [-DEM <str>] [-DIRECTION <str>] [-CONNECTION <str>] [-ORDER <str>] [-BASIN <str>] [-SEGMENTS <str>] [-BASINS <str>] [-NODES <str>] [-THRESHOLD <num>]
## -DEM:<str> Elevation
## Grid (input)
## -DIRECTION:<str> Flow Direction
## Grid (optional output)
## -CONNECTION:<str> Flow Connectivity
## Grid (optional output)
## -ORDER:<str> Strahler Order
## Grid (optional output)
## -BASIN:<str> Drainage Basins
## Grid (optional output)
## -SEGMENTS:<str> Channels
## Shapes (output)
## -BASINS:<str> Drainage Basins
## Shapes (output)
## -NODES:<str> Junctions
## Shapes (optional output)
## -THRESHOLD:<num> Threshold
## Integer
## Minimum: 1rsaga.geoprocessor(lib = 'ta_channels',module = 5,
param = list(DEM = 'Insumos/DEM_fill.sdat',
ORDER = 'Salida/Strahler_order',
SEGMENTS = 'Salida/red_drenaje'),
env = env)##
##
## SAGA Version: 2.3.1
##
## library path: /usr/lib/x86_64-linux-gnu/saga/
## library name: libta_channels
## library : Channels
## tool : Channel Network and Drainage Basins
## author : O.Conrad (c) 2003
## processors : 4 [4]
##
## loading: DEM_fill
##
##
## Parameters
##
##
## Grid system: 30; 558x 661y; 963390.8124x 8058312.3304y
## Elevation: DEM_fill
## Flow Direction: <not set>
## Flow Connectivity: <not set>
## Strahler Order: Strahler Order
## Drainage Basins: <not set>
## Channels: Channels
## Drainage Basins: Drainage Basins
## Junctions: <not set>
## Threshold: 5
##
## Flow Direction
## Stream Order
## Junctions
## Drainage Basins
## Vectorising Grid Classes
##
##
## Parameters
##
##
## Grid system: 30; 558x 661y; 963390.8124x 8058312.3304y
## Grid: Drainage Basins
## Polygons: Polygons
## Class Selection: all classes
## Vectorised class as...: one single (multi-)polygon object
## Keep Vertices on Straight Lines: no
##
## class identification
## Create index: Drainage Basins
## edge detection
## edge collection
## Channels
## Save grid: Salida/Strahler_order...
## Save shapes: Salida/red_drenaje...Finalmente visualizamos nuestros resultados
** Strahler order
raster('Salida/Strahler_order.sdat') %>%
tm_shape() +
tm_basemap('Esri.WorldTopoMap') +
tm_raster(style = 'fixed',breaks = seq(0,5,1),palette = '-viridis')Red de drenaje
read_sf('Salida/red_drenaje.shp') %>%
tm_shape() +
tm_basemap('Esri.WorldTopoMap') +
tm_lines()